/*
** SymmetricLevinson.cs
**
** Copyright (c) 2003-2006, dnAnalytics. All rights reserved.
*/

#region Using directives

using System;
using dnAnalytics.Math;
using dnAnalytics.Exceptions;
using dnAnalytics.Resources;

#endregion Using directives

namespace dnAnalytics.LinearAlgebra.Decomposition
{

    /// <summary>
    /// A Levinson decomposition for symmetric square Toeplitz systems of <b>double</b> type.
    /// </summary>
    /// <remarks>
    /// This class provides members for inverting a symmetric square Toeplitz matrix
    /// (see <see cref="Inverse()"/> member) and calculating the determinant of the matrix
    /// (see <see cref="Determinant()"/> member).
    /// <para>
    /// The class implements an <b>UDL</b> decomposition of the inverse of the Toeplitz matrix.
    /// The decomposition is based upon Levinson's algorithm. As a consequence, all operations
    /// require approximately <b>N</b> squared FLOPS, where <b>N</b> is the matrix order. This
    /// is significantly faster than Cholesky's factorization for symmetric matrices, which
    /// requires <b>N</b> cubed FLOPS.
    /// </para>
    /// <para>
    /// A requirement of Levinson's algorithm is that all the leading sub-matrices and the principal
    /// matrix must be non-singular. During the decomposition, sub-matrices and the principal matrix are
    /// checked to ensure that they are non-singular. When a singular matrix is found, the decomposition
    /// is halted and an internal flag is set. The <see cref="IsSingular"/> property may be used to access
    /// the flag, to determine if any singular matrices were detected.
    /// </para>
    /// <para>
    /// It has been shown that Levinson's algorithm is weakly numerically stable for positive-definite
    /// Toeplitz matrices. It is usual to restrict the use of the algorithm to such
    /// matrix types. The <see cref="IsPositiveDefinite"/> property may be checked to verify that the
    /// Toeplitz matrix is positive-definite.
    /// </para>
    /// <para>
    /// If one of the leading sub-matrices or the principal matrix is near-singular, then the accuracy
    /// of the decomposition will be degraded. An estimate of the resulting error is provided and may be
    /// accessed with the <see cref="Error"/> property. This estimate is only valid for positive-definite
    /// matrices.
    /// </para>
    /// <para>
    /// A outline of this approach to the UDL decomposition of inverse Toeplitz matrices is
    /// found in the following reference:
    /// </para>
    /// <para>
    /// <i>Sun-Yuan Kung and Yu Hen Hu</i>, A Highly Concurrent Algorithm and Pipelined
    /// Architecture for Solving Toeplitz Systems, IEEE Transactions on Acoustics,
    /// Speech and Signal Processing, Volume ASSP-31, Number 1, February 1983, pages
    /// 66 - 75.
    /// </para>
    /// </remarks>
    public class SymmetricLevinson : Algorithm
    {
        #region Fields
        double[] leftColumn;
        double[] diagonal;
        double[][] lowerTriangle;
        int order;
        bool isSingular;
        bool isPositiveDefinite;
        private double det = double.MinValue;

        /// <summary>
        /// The Cybenko error-bound
        /// </summary>
        private double cybenkoError;

        #endregion Fields

        #region Constructors

        /// <summary>
        /// Constructor with <b>Vector</b> parameter.
        /// </summary>
        /// <param name="T">The left-most column of the Toeplitz matrix.</param>
        /// <exception cref="ArgumentNullException">
        /// <b>T</b> is a null reference.
        /// </exception>
        /// <exception cref="RankException">
        /// The length of <b>T</b> is zero.
        /// </exception>
        public SymmetricLevinson(Vector T)
        {
            // check parameter
            if (T == null)
            {
                throw new System.ArgumentNullException("T", Strings.NullParameterException);
            }
            else if (T.Count == 0)
            {
                throw new RankException(Strings.LengthOfZero);
            }

            order = T.Count;
            leftColumn = T.ToArray();

            // allocate memory for lower triangular matrix
            lowerTriangle = new double[order][];
            for (int i = 0; i < order; i++)
            {
                lowerTriangle[i] = new double[i + 1];
            }

            // allocate memory for diagonal
            diagonal = new double[order];
        }

        /// <summary>
        /// Constructor with <b>double</b> array parameter.
        /// </summary>
        /// <param name="T">The left-most column of the Toeplitz matrix.</param>
        /// <exception cref="ArgumentNullException">
        /// <b>T</b> is a null reference.
        /// </exception>
        /// <exception cref="RankException">
        /// The length of <b>T</b> is zero.
        /// </exception>
        public SymmetricLevinson(double[] T)
        {
            // check parameter
            if (T == null)
            {
                throw new System.ArgumentNullException("T", Strings.NullParameterException);
            }
            else if (T.Length == 0)
            {
                throw new RankException(Strings.LengthOfZero);
            }

            order = T.Length;
            leftColumn = new double[order];
            Buffer.BlockCopy(T, 0, leftColumn, 0, order * Constants.SizeOfDouble);

            // allocate memory for lower triangular matrix
            lowerTriangle = new double[order][];
            for (int i = 0; i < order; i++)
            {
                lowerTriangle[i] = new double[i + 1];
            }

            // allocate memory for diagonal
            diagonal = new double[order];
        }
        #endregion Constructors

        #region Properties
        /// <summary>
        /// Get the order of the Toeplitz matrix.
        /// </summary>
        public int Order
        {
            get
            {
                return order;
            }
        }

        #endregion Properties

        #region Public Methods
        /// <summary>
        /// Check if the Toeplitz matrix or any leading sub-matrices are singular.
        /// </summary>
        /// <returns><b>true</b> if the Toeplitz matrix or any leading sub-matrices are singular; 
        /// otherwise <b>false</b>. </returns>
        /// <remarks>
        /// If the Toeplitz matrix or any leading sub-matrices are singular, it is
        /// not possible to complete the <b>UDL</b> decomposition using the Levinson
        /// algorithm.
        /// </remarks>
        public bool IsSingular()
        {
            Compute();
            return isSingular;
        }

        /// <summary>
        /// Check if the Toeplitz matrix is positive definite.
        /// </summary>
        /// <returns><b>true</b> if the Toeplitz matrix is positive definite; otherwise <b>false</b>. </returns>
        /// <remarks>
        /// It has only been shown that the Levinson algorithm is weakly numerically stable
        /// for symmetric positive-definite Toeplitz matrices.
        /// Based on empirical results, it appears that the Levinson algorithm
        /// gives reasonable accuracy for many symmetric indefinite matrices.
        /// It may be desirable to restrict the use of this class to positive-definite matrices
        /// to guarantee accuracy.
        /// </remarks>
        public bool IsPositiveDefinite()
        {
            Compute();
            return isPositiveDefinite;

        }

        /// <summary>
        /// Get an error estimate for the inverse matrix.
        /// </summary>
        /// <returns> An error estimate for the inverse matrix.</returns>
        /// <remarks>
        /// This estimate is approximate and is only valid for positive definite matrices.
        /// It is useful for checking the accuracy of the <b>UDL</b> decomposition. If the
        /// Toeplitz matrix or one of the leading sub-matrices is near-singular,
        /// the accuracy of the decomposition will be degraded. This error estimate
        /// can be used to identify such situations.
        /// <para>
        /// The estimate is based on the Cybenko error-bound for the inverse
        /// and is a relative error. When calculating this error it is assumed that 
        /// the <see cref="double"/> type utilises a 52-bit mantissa.
        /// </para>
        /// <para>
        /// The error estimate will range in value between 2.22e-16 and <see cref="double.PositiveInfinity"/>.
        /// </para>
        /// </remarks>
        public double Error()
        {
            Compute();
            return cybenkoError;
        }

        ///<summary>
        /// Get the lower triangle matrix of the UDL factorisation.
        /// </summary>
        /// <returns>The lower triangle matrix of the UDL factorisation.</returns>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the lower triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public Matrix Lower()
        {
            Compute();
            if (isSingular)
            {
                throw new SingularMatrixException();
            }
            Matrix ret = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Lower(ret);
            return ret;
        }

        ///<summary>
        /// Get the lower triangle matrix of the UDL factorisation.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the lower triangle matrix 
        /// of the UDL factorisation into.</param>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the lower triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void Lower(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            for (int i = 0; i < order; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    result.ValueAt(i, j, lowerTriangle[i][j]);
                }
            }
        }

        ///<summary>
        /// Get the diagonal matrix of the UDL factorisation.</summary>
        /// <remarks>
        /// <returns>The diagonal matrix of the UDL factorisation.</returns>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the diagonal matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public Matrix Diagonal()
        {
            Compute();
            if (isSingular)
            {
                throw new SingularMatrixException();
            }
            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Diagonal(result);
            return result;
        }

        ///<summary>
        /// Get the diagonal matrix of the UDL factorisation.</summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the diagonal matrix 
        /// of the UDL factorisation into.</param>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the diagonal matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void Diagonal(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            for (int i = 0; i < diagonal.Length; i++)
            {
                result.ValueAt(i, i, diagonal[i]);
            }

        }

        ///<summary>
        /// Get the upper triangle matrix of the UDL factorisation.</summary>
        /// <remarks>
        /// <returns>The upper triangle matrix of the UDL factorisation.</returns>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the upper triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public Matrix Upper()
        {
            Compute();
            if (isSingular)
            {
                throw new SingularMatrixException();
            }
            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Upper(result);
            return result;
        }

        /// <summary>
        /// Get the upper triangle matrix of the UDL factorisation.</summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the upper triangle 
        /// matrix of the UDL factorisation into.</param>
        /// <remarks>It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the upper triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void Upper(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            for (int i = 0; i < order; i++)
            {
                for (int j = i; j < order; j++)
                {
                    result.ValueAt(i, j, lowerTriangle[j][i]);
                }
            }
        }

        /// <summary>
        /// Get a copy of the Toeplitz matrix.
        /// </summary>
        /// <returns>A copy of the Toeplitz matrix.</returns>
        public Matrix ToeplitzMatrix()
        {
            Compute();
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            ToeplitzMatrix(result);
            return result;
        }

        /// <summary>
        /// Get a copy of the Toeplitz matrix.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the Toeplitz matrix into.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void ToeplitzMatrix(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            if (order > 1)
            {
                for (int j = 0; j < leftColumn.Length; j++)
                {
                    result.ValueAt(0, j, leftColumn[j]);
                }

                // fill bottom row (reverse order)
                double[] bottom = new double[order];

                for (int i = 0, j = order - 1; i < order; i++, j--)
                {
                    bottom[i] = leftColumn[j];
                }

                double[] temp = new double[order];

                // fill rows in-between
                for (int i = 1, j = order - 1; j > 0; i++)
                {
                    Buffer.BlockCopy(leftColumn, 0, temp, i * Constants.SizeOfDouble, (j--) * Constants.SizeOfDouble);
                    Buffer.BlockCopy(bottom, j * Constants.SizeOfDouble, temp, 0, i * Constants.SizeOfDouble);
                    result.SetRow(i, temp);
                }
            }
            else
            {
                result[0, 0] = leftColumn[0];
            }
        }

        /// <summary>
        /// Get the determinant of the Toeplitz matrix.
        /// </summary>
        /// <returns>The determinant</returns>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the determinant.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public double Determinant()
        {
            Compute();
            // if any of the matrices is singular give up
            if (isSingular)
            {
                return 0;
            }
            else
            {
                if (det == double.MinValue)
                {
                    det = 1.0;
                    for (int i = 0; i < order; i++)
                    {
                        det *= diagonal[i];
                    }
                    det = 1.0 / det;
                }
            }

            return det;
        }

        /// <summary>
        /// Get the inverse of the Toeplitz matrix.
        /// </summary>
        /// <returns>The inverse matrix.</returns>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before using an extended version of Trench's algorithm to complete
        /// the inversion.
        /// <para>
        /// The extended version of Trench's algorithm requires approximately <b>N</b> squared FLOPS,
        /// compared to <b>N</b> cubed FLOPS if we simply multiplied the <b>UDL</b> factors
        /// (<b>N</b> is the matrix order).
        /// </para>
        /// </remarks>
        public Matrix Inverse()
        {
            Compute();
            if (isSingular)
            {
                throw new SingularMatrixException();
            }
            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Inverse(result);
            return result;
        }

        /// <summary>
        /// Get the inverse of the Toeplitz matrix.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the inverse of the Toeplitz matrix into.</param>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        /// <remarks>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before using an extended version of Trench's algorithm to complete
        /// the inversion.
        /// <para>
        /// The extended version of Trench's algorithm requires approximately <b>N</b> squared FLOPS,
        /// compared to <b>N</b> cubed FLOPS if we simply multiplied the <b>UDL</b> factors
        /// (<b>N</b> is the matrix order).
        /// </para>
        /// </remarks>
        public void Inverse(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            double[] A = lowerTriangle[order - 1];
            double A1, A2, scale;
            int i, j, k, l;

            // setup the first row in wedge
            scale = diagonal[order - 1];
            for (i = 0, j = order - 1; i < order; i++, j--)
            {
                result.ValueAt(0, i, scale * A[j]);
            }

            // calculate values in the rest of the wedge
            for (i = 1; i < (1 + order) / 2; i++)
            {
                A1 = A[order - i - 1];
                A2 = A[i - 1];
                for (j = i, k = i - 1, l = order - i - 1; j < order - i; j++, k++, l--)
                {
                    double temp = result.ValueAt(i - 1, k) + scale * (A1 * A[l] - A2 * A[k]);
                    result.ValueAt(i, j, temp);
                }
            }
            // this is symmetric matrix ...
            for (i = 0; i < (1 + order) / 2; i++)
            {
                for (j = i; j < order - i; j++)
                {
                    result.ValueAt(j, i, result.ValueAt(i, j));
                }
            }

            // and a persymmetric matrix.
            for (i = 0, j = order - 1; i < order; i++, j--)
            {
                for (k = 0, l = order - 1; k < j; k++, l--)
                {
                    result.ValueAt(l, j, result.ValueAt(i, k));
                }
            }
        }

        /// <summary>
        /// Solve a symmetric square Toeplitz system with a right-side <see cref="Vector"/>.
        /// </summary>
        /// <param name="Y">The right-hand side of the system.</param>
        /// <returns>The solution vector.</returns>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="Y"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// The length of <paramref name="Y"/>  is not equal to the number of rows in the Toeplitz matrix.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// the symmetric square Toeplitz matrix, <b>X</b> is the unknown solution vector
        /// and <b>Y</b> is a known vector.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, and then calculates the solution vector.
        /// </para>
        /// </remarks>
        public Vector Solve(Vector Y)
        {
            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }
            if (order != Y.Count)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }
            Vector X = VectorBuilder.CreateVector(order, VectorType.Dense);
            Solve(Y, X);
            return X;
        }
   
        /// <summary>
        /// Solve a symmetric square Toeplitz system with a right-side <see cref="Vector"/>.
        /// </summary>
        /// <param name="Y">The right-hand side of the system.</param>
        /// <param name="result">The solution vector.</param>
        /// <exception cref="ArgumentNullException">
        /// If either <paramref name="Y"/> or <paramref name="result"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// If the length of <paramref name="Y"/> or <paramref name="result"/> is not equal to the number of rows in the Toeplitz matrix.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// the symmetric square Toeplitz matrix, <b>X</b> is the unknown solution vector
        /// and <b>Y</b> is a known vector.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, and then calculates the solution vector.
        /// </para>
        /// </remarks>
        public void Solve(Vector Y, Vector result)
        {
            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }
            if (order != Y.Count)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            if (result == null)
            {
                throw new System.ArgumentNullException("result", Strings.NullParameterException);
            }
            if (order != result.Count)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            
            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            int i, j, l;			// index/loop variables
            double Inner;			// inner product
            double G;				// scaling constant
            double[] A;				// reference to current order coefficients

            // setup zero order solution
            result[0] = Y[0] / leftColumn[0];

            // solve systems of increasing order
            for (i = 1; i < order; i++)
            {
                // calculate inner product
                Inner = Y[i];
                for (j = 0, l = i; j < i; j++, l--)
                {
                    Inner -= result[j] * leftColumn[l];
                }
                // get the current predictor coefficients row
                A = lowerTriangle[i];

                // update the solution vector
                G = Inner * diagonal[i];
                for (j = 0; j <= i; j++)
                {
                    result[j] += G * A[j];
                }
            }
        }

        /// <summary>
        /// Solve a symmetric square Toeplitz system with a right-side matrix.
        /// </summary>
        /// <param name="Y">The right-hand side of the system.</param>
        /// <returns>The solution matrix.</returns>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="Y"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// The number of rows in <paramref name="Y"/> is not equal to the number of rows in the Toeplitz matrix.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// a symmetric square Toeplitz matrix, <b>X</b> is the unknown solution matrix
        /// and <b>Y</b> is a known matrix.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, and then calculates the solution matrix.
        /// </para>
        /// </remarks>
        public Matrix Solve(Matrix Y)
        {
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }
            if (order != Y.Rows)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }
            DenseMatrix X = (DenseMatrix)MatrixBuilder.CreateMatrix(order, Y.Rows, MatrixType.Dense);
            Solve(Y, X);
            return X;
        }

        /// <summary>
        /// Solve a symmetric square Toeplitz system with a right-side matrix.
        /// </summary>
        /// <param name="Y">The right-hand side of the system.</param>
        /// <param name="result">The solution matrix.</param>
        /// <exception cref="ArgumentNullException">
        /// If either <paramref name="Y"/> or <paramref name="result"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// <para>The number of rows in <paramref name="Y"/> is not equal to the number of rows in the Toeplitz matrix.</para>
        /// <para>- or -</para>
        /// <para>The number of rows in <paramref name="result"/> is not equal to the number of columns in the Toeplitz matrix.</para>
        /// <para>- or -</para>
        /// <para>The number of rows in <paramref name="Y"/> is not equal to the number of columns in <paramref name="result"/>.</para>
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// a symmetric square Toeplitz matrix, <b>X</b> is the unknown solution matrix
        /// and <b>Y</b> is a known matrix.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, and then calculates the solution matrix.
        /// </para>
        /// </remarks>        
        public void Solve(Matrix Y, Matrix result)
        {
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }
            if (order != Y.Rows)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }
            
            if (result == null)
            {
                throw new System.ArgumentNullException("result", Strings.NullParameterException);
            }
            if (order != result.Rows)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }
            if (Y.Columns != result.Columns)
            {
                throw new NotConformableException(Strings.ParametersNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            int M = Y.Rows;
            int i, j, l, m;			// index/loop variables
            double[] Inner;			// inner product
            double[] G;				// scaling constant
            double[] A;				// reference to current order coefficients
            double scalar;
            double temp;

            // allocate memory for solution

            Inner = new double[M];
            G = new double[M];

            // setup zero order solution
            scalar = 1.0 / leftColumn[0];
            for (m = 0; m < M; m++)
            {
                result.ValueAt(0, m, (scalar * Y.ValueAt(0, m)));
            }

            // solve systems of increasing order
            for (i = 1; i < order; i++)
            {
                // calculate inner product
                for (m = 0; m < M; m++)
                {
                    Inner[m] = Y[i, m];
                }

                for (j = 0, l = i; j < i; j++, l--)
                {
                    scalar = leftColumn[l];
                    for (m = 0; m < M; m++)
                    {
                        Inner[m] -= scalar * result.ValueAt(j, m);
                    }
                }

                // get the current predictor coefficients row
                A = lowerTriangle[i];

                // update the solution matrix
                for (m = 0; m < M; m++)
                {
                    G[m] = diagonal[i] * Inner[m];
                }
                for (j = 0; j <= i; j++)
                {
                    scalar = A[j];
                    for (m = 0; m < M; m++)
                    {
                        temp = result.ValueAt(j, m) + scalar * G[m];
                        result.ValueAt(j, m, temp);
                    }
                }
            }
        }
        #endregion Public Methods

        #region Protected Methods
        /// <summary>
        /// Calculate the UDL decomposition of the inverse Toeplitz matrix.
        /// </summary>
        /// <remarks>
        /// The member checks each iteration for a singular principal sub-matrix.
        /// If a singular sub-matrix is detected, the private boolean <i>m_IsSingular</i>
        /// is set to <b>true</b>, the <i>m_IsPositiveDefinite</i>EM> flag is set <b>false</b> and the
        /// decomposition halted.
        /// <para>
        /// When the decomposition is completed, the lower triangular matrix is stored
        /// in the jagged array <i>m_LowerTriangle</i>, the diagonal elements
        /// in the array <i>m_Diagonal</i> and the private boolean <i>m_IsSingular</i>
        /// is set to <b>false</b>.
        /// </para>
        /// <para>
        /// The diagonal elements are checked to determine if the Toeplitz matrix is
        /// positive definite. The <i>m_IsPositiveDefinite flag</i>EM> is set <b>true</b>, if all
        /// diagonal elements are positive, otherwise it is set to <b>false</b>.
        /// </para>
        /// <para>
        /// During the calculation the reflection coefficients are used to calculate the
        /// Cybenko error estimate which is stored in the member <i>m_CybenkoError</i>.
        /// </para>
        /// </remarks>
        protected override void InternalCompute()
        {
            int i, j, l;			// index/loop variables
            double Inner;			// inner product
            double K;				// reflection coefficient
            double[] B;				// reference to previous order coefficients
            double[] A;				// reference to current order coefficients


            // check if principal diagonal is zero
            if (leftColumn[0] == 0.0)
            {
                isSingular = true;
                isPositiveDefinite = false;
                cybenkoError = double.MaxValue;
                return;
            }

            // setup zero order solution
            B = lowerTriangle[0];
            B[0] = 1.0;
            diagonal[0] = 1.0 / leftColumn[0];
            cybenkoError = Constants.UnitRoundOff;

            // solve systems of increasing order
            for (i = 1; i < order; i++, B = A)
            {
                // calculate inner product
                Inner = 0.0;
                for (j = 0, l = 1; j < i; j++, l++)
                {
                    Inner += leftColumn[l] * B[j];
                }

                // calculate the reflection coefficient
                K = -Inner * diagonal[i - 1];

                // get the current low triangle row
                A = lowerTriangle[i];

                // update low triangle elements
                A[0] = 0.0;
                Buffer.BlockCopy(B, 0, A, 1 * Constants.SizeOfDouble, i * Constants.SizeOfDouble);
                //Array.Copy(B, 0, A, 1, i);
                for (j = 0, l = i - 1; j < i; j++, l--)
                {
                    A[j] += K * B[l];
                }

                // check if singular sub-matrix
                if (K * K == 1.0)
                {
                    isSingular = true;
                    isPositiveDefinite = false;
                    cybenkoError = double.MaxValue;
                    return;
                }

                // update diagonal
                diagonal[i] = diagonal[i - 1] / (1.0 - K * K);

                // update error estimate
                cybenkoError *= (1.0 + System.Math.Abs(K)) / (1.0 - System.Math.Abs(K));
            }

            // check if system is positive-definite
            for (i = 0, isPositiveDefinite = true; (i < order) && (isPositiveDefinite); i++)
            {
                if (diagonal[i] < 0.0)
                {
                    isPositiveDefinite = false;
                }
            }

            // solution completed
            isSingular = false;
            return;
        }

        #endregion Protected Methods
    }
}
